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Abstract 

Little attention has been paid to the use of multi-sample batch-marking studies, 
as it is generally assumed that an individual's capture history is necessary for 
fully efficient estimates. However, recently, Huggins et al. (2010) present a 
pseudo-likelihood for a multi-sample batch-marking study where they used 
estimating equations to solve for survival and capture probabilities and then 
derived abundance estimates using a Horvitz-Thompson-type estimator. We 
have developed and maximized the likelihood for batch-marking studies. We 
use data simulated from a Jolly-Seber-type study and convert this to what 
would have been obtained from an extended batch-marking study. We compare 
our abundance estimates obtained from the Crosbie-Manly-Arnason-Schwarz 
(CMAS) model with those of the extended batch-marking model to determine 
the efficiency of collecting and analyzing batch-marking data. We found that 
estimates of abundance were similar for all three estimators: CMAS, Huggins, 
and our likelihood. Gains are made when using unique identifiers and employ- 
ing the CMAS model in terms of precision; however, the likelihood typically 
had lower mean square error than the pseudo-likelihood method of Huggins 
et al. (2010). When faced with designing a batch-marking study, researchers 
can be confident in obtaining unbiased abundance estimators. Furthermore, 
they can design studies in order to reduce mean square error by manipulating 
capture probabilities and sample size. 



Introduction 

Batch-marking experiments have largely been neglected 
by statistical ecologists, as they are deemed inferior and 
to be avoided (Pollock 1981; Pollock and Mann 1983). 
However, biologists still use batch-marking for various 
purposes, and for some studies, they may be the only 
option available (e.g., insects, juvenile fish). 

There are other types of batch-marking studies that dif- 
fer in design from the one we study here. For example, 
Measey et al. (2003) performed a three-sample batch- 
marking experiment on caecilians where individuals were 
given a batch mark on the first occasion and on the sec- 
ond occasion, a subsample was given a secondary mark. 
Both marked and unmarked captured individuals were 
recorded at each sample time. Because an individual's 
capture history can be deduced when a different batch 
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mark is applied on each sampling occasion, a Jolly-Seber 
type model can be fitted to analyze these data (Jolly 1965; 
Seber 1965). However, the disadvantage of the design is 
that there is a physical limitation to how many marks can 
be applied to an individual and this would vary by both 
species and mark type. 

Frequently, batch marks are used to study movement 
of individuals between locations. For example, Roberts 
and Angermeier (2007) studied the movements of three 
fish species in the South Fork of the Roanoke River, Vir- 
ginia, using a two-sample study. Here, they constructed 
movement corridors with favorable pool characteristics 
between suitable habitats and compared movement rates 
in corridors with unfavorable characteristics. Captured 
fish were given a mark that was a randomly assigned 
color and body location. Recaptured individuals were 
counted, and movement rates were estimated. 
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Skalski et al. (2009) review several batch-marking 
designs and marking methods for very small fish. How- 
ever, most of these result in complete capture history 
information and are for two or three sampling occasion 
designs. 

Arguments against using batch marks are based on the 
lack of individual capture histories. For example, if a 
marked individual is captured at sample time three, it is 
not known whether this individual was one of the marked 
individuals captured at sample time two or not. In addi- 
tion, batch-marking experiments do not allow for ade- 
quate testing of model assumptions (Williams et al. 2002; 
p. 312). 

We motivate this work with the data found in Huggins 
et al. (2010). They describe a study of oriental weatherlo- 
ach (Misgurnus anguillicaudatu), which is a freshwater fish 
native to Eurasia and Northern Africa. It was brought to 
Australia for use in aquaria but was accidentally released, 
and the aim of the study was to investigate activity pat- 
terns of the wild populations in the inland waters. 

Huggins et al. (2010) provide a pseudo-likelihood 
method for analyzing an extended batch-marking study. 
They caution that the likelihood is intractable as the 
number of marked individuals alive at any sample time 
is unknown. They condition on released individuals to 
develop estimating equations and obtain capture and 
survival probability estimates. Then, they use a Horvitz- 
Thompson-type estimator to estimate population size at 
each time point after obtaining capture probability esti- 
mates. Standard errors are obtained by first using a 
sandwich estimator for the variance of the model 
parameters (Freedman 2012) and then using the delta 
method to obtain estimated standard errors for popula- 
tion size. 

We develop the batch-marking likelihood conditional 
on release (rather than the pseudo-likelihood), followed 
by a Horvitz-Thompson-like estimator for abundance. 
Although theoretically the likelihood can be maximized, it 
involves nested summations, resulting in a large number 
of computations, but the calculations can be run in paral- 
lel when a multiprocessor computer, a cluster or a grid is 
available. For this article, we investigate the use of 
extended batch-marking data in comparison with the 
Crosbie-Manly-Arnason-Schwarz (CMAS) model (Schwarz 
and Arnason 1996) to study the loss in estimation preci- 
sion when one does not have information on individual 
encounter histories for a seven sampling occasion 
simulation experiment under various parameter values. 

Materials and Methods 

An extended batch-marking study is one where individu- 
als captured at the first sample time are all given the same 



nonunique type of tag (e.g., blue in color). At subsequent 
sample times, individuals captured with tags are counted 
and unmarked individuals are given a different color 
batch mark resulting in an independent cohort. Table 1 
provides an example of generated data from a four sam- 
pling occasion extended batch-marking experiment. New 
marks are not given to marked individuals, and thus indi- 
vidual capture histories cannot be obtained. For example, 
it is not known whether the nine blue-tagged individuals 
at sample time three are a subset of the 16 found at time 
two. Note the similarity to the m-array notation for 
Cormack-Jolly-Seber data (see Williams et al. 2002; 
p. 419). 

The assumptions we make are similar to other open 
population capture-recapture models namely: 

• All individuals behave independently. 

• All individuals have the same probability of capture at 
sample time j; j = 1, 2, k. 

• All individuals have the same probability of survival 
between sample times; and j + 1; j = 1, 2, k — 1. 

• Individuals do not lose their tags. 

Below we detail notation used in the model develop- 
ment. 

Statistics or indices 

i index for release occasion (or colour of tag). 

j index for recapture occasion. 

k the number of sampling occasions. 

Tij the number of individuals tagged and released at time 

x and recaptured at time j, i = 1, 2, k — 1; 

j = i + 1, k. 
Ri the number of individuals released at time ;'; i = 1, 2, 
k. 

Latent variables 

My the number of marked individuals released at sample 
time i, alive and available for capture at sample time 
j; j = i k. Note that M u = R,. 



Table 1. Example data for the extended bmarking design. The num- 
ber of individuals marked with a particular tag color is found on the 
diagonal, while the number of recaptures is on the off diagonal. 



Release 


Color 


Occasion 






1 2 


3 


4 


1 


Blue 


21 16 


9 


11 


2 


Green 


22 


15 


12 


3 


Orange 




17 


4 


4 


Red 






6 
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djj the number of deaths between sample times j and j 
+ 1 from release group z; £?,)' = Mu - Mj,- +J ; 
i — 1, k; j = i, k — 1. 

Parameters 

0y the probability of survival for individuals from release 
group i between times ; and ; + l;j = 1, 2, k— 1. 
the probability of capture for individuals from release 
group z' at time j; j = 2, 3, k. 

We develop the likelihood by first looking at the joint 
distribution of the recaptures r i; - and deaths d y - given the 
releases Rj. We then obtain the marginal distribution of 
the recaptures given releases by summing over all possible 
values of the deaths. The likelihood can be written as 

k k 

i (^^) = IlE---E t[ p to\*»,-.;di*-uBi) 

1=1 d„ lfcj_ij=t+l (1) 

P(du,...,di,k-i\Ri)- 

Conditional on release, we model the recaptures as 
independent given deaths r-,j \ dn , . . . , d lk _\ , -R; ~ Binomial 
(Ri ~Ylm=i d im,Pij) and the deaths as d ib ...,d ik \Ri~ 
Multinomial(Rj, n ih ... ,n ik ) where ^ = (1 - Ilm=i 
0 im . We note that 71* = 1 - Ylmli n im and d ik = R~ 
Y^m=i dim where d ik would be the individuals that were 
released at time i and are still alive after the last sample 
time k. These d ik are convenient for modeling purposes. 
Thus, the likelihood becomes 

i=l d,j 




Inference 

The calculation of the likelihood involves nested summa- 
tions for the latent d^ variables which require high execution 
times if serially computed or cause the available RAM to be 
used up if fully vectorized. We developed parallel computer 
code to implement this model in MATLAB (MATLAB 
2012) trading off CPU speed and memory that works for up 
to 1 1 sampling occasions. For experiments beyond 1 1 sam- 
pling occasions, we propose to use our likelihood up to the 
1 1th sample time; and the pseudo-likelihood (Huggins et al. 



2010) for occasions 12 through k. In the simulation studies, 
we first produce maximum likelihood estimates of the sur- 
vival and capture probability parameters. Then, we derive a 
Horvitz-Thompson-type estimator for population size at 
each sample time (Nj) using the capture probability esti- 
mates and the number of individuals captured at each sam- 
pling time, Nj = J2i r ij/Pj- Standard errors for the Nj are 
estimated from the estimated variances of r y - and pj using 
the delta method (see Huggins et al. 2010 for details). 

Monte Carlo simulations 

We simulated data from a k = 7 sample occasion Cros- 
bie-Manly-Arnason-Schwarz model (Schwarz and Arna- 
son 1996) with constant survival probabilities (0 = 0.2, 
0.5, or 0.8), constant capture probabilities (p = 0.2, 0.5, 
or 0.8), and entry probabilities equal across time (l/k) 
with both a small superpopulation size (N = 200) and a 
larger superpopulation size (N = 1000). The superpopula- 
tion N is defined as the population that enters the popu- 
lation of interest at some point during the study. The 
computing time (on a dual quad core 2.53 GHz, 32 Gb 
RAM Linux server) was approximately 2 days to run 100 
replications for N = 1000 but larger superpopulation size, 
more samples, or more replications would result in longer 
computing times. Parameter values were selected to 
obtain sparse-to-plentiful data by varying the probability 
of capture and survival. 

For each set of parameter values, we simulated 100 
CMAS datasets and collapsed these datasets into batch- 
marking data. We analyzed the individual capture history 
data using the CMAS model with constant parameters 
implemented in RMark (Laake 2013). The associated 
batch-marking data were analyzed using both the 
pseudo-likelihood (Huggins et al. 2010) and the likeli- 
hood with constant parameters (4>, p). For all methods of 
analyses, we estimated the survival and capture probabil- 
ities and obtained abundance estimates and estimated 
standard errors for each sampling time. The 100 dataset 
results are summarized using box plots for the estimated 
abundance, and estimated capture and survival probabili- 
ties. Root mean square errors (v/ 100 -1 Yll=i@r ~ @) 2 ) 
were calculated for the capture and survival probabilities. 
While standard errors were estimated, plots of these results 
are not included in the interest of space but are provided 
in the supplementary materials (see Supporting 
information). 

Results 

Figures 1 and 2 provide results for the 9 simulation stud- 
ies under varying parameter values for estimates of p and 
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Fig. 1. Boxplots of capture probability 
estimates (p) from 100 simulated datasets, for 
the Crosbie-Manly-Arnason-Schwarz (C: red), 
the likelihood (L: green), and the pseudo- 
likelihood (H: yellow; Huggins et al. 2010) 
methods when parameter values are W = 200, 
and p = 0.2, 0.5, 0.8 for <j> = 0.2 (top), 
(j> = 0.5 (middle), and </> = 0.8 (bottom). 




0.8 P = 0.2 



P= 0.5 



P = 0.8 



cp, respectively, for N = 200. For sparse data (p = 0.2 or 
p = 0.5, cp = 0.2), many of the simulations produced 
boundary estimates for p and occasionally cp (see Fig. 1), 
and the calculation of the standard errors then failed due 
to the Hessian being singular. These simulation failures 
are similar to what happens in the Cormack-Jolly-Seber 
model when analytical estimates of <f> exceed 1 with sparse 
data and actual parameter values are close to 0 and 1. In 
these cases, the maximization function in MATLAB con- 
strained estimates to be admissible, i.e. between 0 and 1 
(inclusive). When the estimation of p was on a boundary, 



Nj was estimated at infinity (p = 0) or Y,i r ij (P = !)• 
Table 2 provides the number of simulations out of 100 
that produced boundary estimates for p or cp for all three 
methods. Similar figures for N = 1000 are provided in the 
Supporting information. Results for the estimation of 
standard errors are based on those simulations that did 
not fail (see Supporting information). 

The root mean square error (RMSE) for estimates of p 
and <f> are given in Tables 3 and 4, respectively. As 
expected, we find that within a method, RMSE decreases 
as p and (p increase. Similarly, RMSE decreases with 
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Fig. 2. Boxplots of survival probability 
estimates (<j>) from 100 simulated datasets, for 
the Crosbie-Manly-Arnason-Schwarz (C: red), 
the likelihood (L: green), and the pseudo- 
likelihood (H: yellow; Huggins et al. 2010) 
methods when parameter values are N = 200, 
and 4> = 0.2, 0.5, 0.8 for p = 0.2 (top), 
p = 0.5 (middle), and p = 0.8 (bottom). 



increased N. We also confirm that the CMAS method 
typically has lower RMSE than either of the batch mark- 
ing methods and that the likelihood method typically has 
lower RMSE than the pseudo-likelihood method. Excep- 
tions to this occur with sparse data when the estimates 
are not reliable. 

Under sparse data conditions (e.g., N = 200, </> = 0.8, 
p = 0.2), the average population size estimates are similar 
between the three methods; however, variability in 
estimates is higher for the likelihood and pseudo-likelihood 
methods as expected (Fig. 3; box plots for other sets of 



parameters are provided in Supporting information). For 
example, average population size estimates for time three 
were 71, 74, and 72 individuals for the CMAS, likelihood, 
and pseudo-likelihood, respectively, and the corresponding 
average standard error estimates were 18, 32, and 32 indi- 
viduals. For higher quality data (e.g., N = 1000, p = 0.5, 
<p = 0.8), we found similar results. The CMAS model 
produces more precise estimates followed by the likelihood 
(see Supporting information for box plots of estimated 
standard errors). For example, the average population size 
estimate for sample time three was 348, 349, and 349 indi- 
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Table 2. The number of simulations out of 100 that produced 
boundary estimates (0 or 1) for either parameter p or </> for the Cros- 
bie-Manly-Arnason-Schwarz (CMAS), the likelihood (L), and the 
pseudo-likelihood (H; Huggins et al. 2010) methods under the 18 sim- 
ulation scenarios. 









<t> 




N Method 


P 


0.2 


0.5 


0.8 


200 CMAS 


0.2 


95 


9 


2 




0.5 


32 


0 


0 




0.8 


37 


0 


0 


L 


0.2 


75 


9 


6 




0.5 


32 


2 


0 




0.8 


36 


18 


0 


H 


0.2 


76 


10 


6 




0.5 


26 


2 


0 




0.8 


31 


16 


0 


1000 CMAS 


0.2 


32 


0 


0 




0.5 


1 


0 


1 




0.8 


1 


0 


0 


L 


0.2 


27 


0 


0 




0.5 


3 


0 


0 




0.8 


9 


2 


0 


H 


0.2 


24 


0 


0 




0.5 


8 


0 


0 




0.8 


14 


1 


0 



Table 4. Root mean square error for estimates of <j> for the Crosbie- 
Manly-Arnason-Schwarz (CMAS), the likelihood (L), and the pseudo- 
likelihood (H; Huggins et al. 2010) methods under the 18 simulation 
scenarios. 



<l> 



N Method 


P 


0.2 


0.5 


0.8 


200 CMAS 


0.2 


0.224 


0.196 


0.095 




0.5 


0.087 


0.068 


0.040 




0.8 


0.041 


0.036 


0.027 


L 


0.2 


0.221 


0.197 


0.106 




0.5 


0.101 


0.083 


0.058 




0.8 


0.064 


0.053 


0.035 


H 


0.2 


0.215 


0.206 


0.116 




0.5 


0.112 


0.090 


0.060 




0.8 


0.075 


0.060 


0.039 


1000 CMAS 


0.2 


0.121 


0.072 


0.046 




0.5 


0.038 


0.027 


0.014 




0.8 


0.021 


0.016 


0.009 


L 


0.2 


0.116 


0.093 


0.067 




0.5 


0.049 


0.035 


0.021 




0.8 


0.035 


0.025 


0.016 


H 


0.2 


0.116 


0.096 


0.067 




0.5 


0.053 


0.037 


0.023 




0.8 


0.038 


0.027 


0.017 



Table 3. Root mean square error for estimates of p for the Crosbie- 
Manly-Arnason-Schwarz (CMAS), the likelihood (L), and the pseudo- 
likelihood (H; Huggins et al. 2010) methods under the 18 simulation 
scenarios. 



N Method 


P 


0.2 


0.5 


0.8 


200 CMAS 


0.2 


0.704 


0.282 


0.062 




0.5 


0.319 


0.100 


0.046 




0.8 


0.160 


0.065 


0.028 


L 


0.2 


0.616 


0.271 


0.065 




0.5 


0.330 


0.136 


0.077 




0.8 


0.215 


0.123 


0.058 


H 


0.2 


0.616 


0.296 


0.073 




0.5 


0.309 


0.151 


0.079 




0.8 


0.233 


0.125 


0.068 


1000 CMAS 


0.2 


0.470 


0.066 


0.026 




0.5 


0.149 


0.036 


0.019 




0.8 


0.086 


0.027 


0.011 


L 


0.2 


0.448 


0.083 


0.035 




0.5 


0.185 


0.050 


0.031 




0.8 


0.137 


0.057 


0.031 


H 


0.2 


0.428 


0.087 


0.035 




0.5 


0.212 


0.053 


0.032 




0.8 


0.150 


0.061 


0.034 



viduals for the CMAS, likelihood, and pseudo-likelihood 
method, respectively, with corresponding average estimated 
standard errors of 12, 29, and 29 individuals. 



Discussion 

For an extended batch-marking study, using the likelihood 
provides more accurate estimates and lower standard 
errors than using the pseudo-likelihood method of Hug- 
gins et al. (2010). However, the computing power neces- 
sary to calculate the likelihood by summing over all 
possible values of deaths is prohibitive when sampling 
times go beyond k = 11 or if R; is large. In these cases, the 
pseudo-likelihood method is computationally faster and 
provides unbiased estimates with similar precision to the 
likelihood approach. Ultimately, if full-capture histories 
are possible and available, then naturally the CMAS model 
outperforms both batch-marking models. With plentiful 
data (large numbers of recaptures), our model can have a 
relative efficiency of between about 30-40% compared 
with the CMAS model; thus, using the CMAS model has 
obvious gains in precision. 

In many of the plots for Nj, the average population size 
increases over time. This is due to the models allowing 
for births/immigration into the population from the 
superpopulation (N). In these simulations, entry probabil- 
ities were equal across time and summed to one. Thus, 
with high survival rates, population size would naturally 
increase with time. 

Practitioners who are confined to using batch marks 
should design their studies to have large sample sizes and 
high capture rates so as to minimize mean square error. 
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■on 



CLHCLHCLHCLHCLHCLHCLH 



P=0.2 1 




P=0.5 1 




P=0.8 1 



For future work, we will complete the model develop- 
ment by incorporating both tagged and untagged individ- 
uals at each sample time. We will also deal with issues 
such as goodness of fit, model selection, and parameter 
redundancy. With the many latent variables in the 
complete data, this model lends itself well to Bayesian 
methods where a state-space formulation is under devel- 
opment. 

With permanent batch marks, tag loss would not be an 
issue. However, if injectable color tags are used for example, 
tag loss may bias parameter estimates. If it were possible to 
double tag individuals, an extended batch-marking model 



Fig. 3. Boxplots of abundance estimates (A/y) 
for each sample time (k = 7) from 100 
simulated datasets, for the Crosbie-Manly- 
Arnason-Schwarz (C: red), the likelihood (L: 
green), and the pseudo-likelihood (H: yellow; 
Huggins et al. 2010) methods when parameter 
values are N = 200, and <j> = 0.8 for p = 0.2 
(top), p = 0.5 (middle), and p = 0.8 (bottom). 
The long black horizontal lines show the 
expected population size at time /. 

incorporating tag retention rates could be developed using 
methods similar to Cowen and Schwarz (2006). However, 
for those study species where double tagging is not possible 
(e.g., insects), separate experiments to estimate tag reten- 
tion would have to be carried out and this auxiliary infor- 
mation could be used to adjust parameter estimates using 
methods similar to Arnason and Mills (1981). 
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Supporting Information 

Additional Supporting Information may be found in the 
online version of this article: 

Figure SI. Boxplots of capture probability estimates (p) 
of 100 simulated datasets, for the Crosbie-Manley-Arna- 
son-Schwarz (C: red), the likelihood (L: green), and the 
pseudo-likelihood (H: yellow; Huggins et al. 2010) when 
parameters values are N = 1000, and p = 0.2,0.5,0.8 for 
4> = 0.2 (top), and <p = 0.5 (middle), <p = 0.8 (bottom). 
Figure S2. Boxplots of survival probability estimates <p of 
100 simulated datasets, for the Crosbie-Manley-Arnason- 
Schwarz (C: red), the likelihood (L: green), and the 
pseudo-likelihood (H: yellow; Huggins et al. 2010) when 
parameters values are N = 1000, and <j) = 0.2, 0.5, 0.8 for 
p = 0.2 (top), andp = 0.5 (middle), p = 0.8 (bottom). 
Figure S3. Boxplots of abundance estimates Nj for each 
sample time (k = 7) of 100 simulated datasets, for the 
Crosbie-Manley-Arnason-Schwarz (C: red), the likelihood 
(L: green), and the pseudo-likelihood (H: yellow; Huggins 
et al. 2010) when parameters values are N = 200, and 
4> = 0.2 for p = 0.2 (top), and p = 0.5 (middle), p = 0.8 
(bottom). The long black horizontal lines show the 
expected population size at time j. 

Figure S4. Boxplots of abundance estimates Nj for each 
sample time (k = 7) of 100 simulated datasets, for the 
Crosbie-Manley-Arnason-Schwarz (C: red), the likelihood 
(L: green), and the pseudo-likelihood (H: yellow; Huggins 
et al. 2010) when parameters values are N = 200, and 
4> = 0.5 for p = 0.2 (top), and p = 0.5 (middle), p = 0.8 
(bottom). The long black horizontal lines show the 
expected population size at time j. 

Figure S5. Boxplots of abundance estimates Nj for each 
sample time (k = 7) of 100 simulated datasets, for the 
Crosbie-Manley-Arnason-Schwarz (C: red), the likelihood 
(L: green), and the pseudo-likelihood (H: yellow; Huggins 
et al. 2010) when parameters values are N = 1000, and 
4> = 0.2 for p = 0.2 (top), and p = 0.5 (middle), p = 0.8 
(bottom). The long black horizontal lines show the 
expected population size at time j. 

Figure S6. Boxplots of abundance estimates Nj for each 
sample time (k = 7) of 100 simulated datasets, for the 
Crosbie-Manley-Arnason-Schwarz (C: red), the likelihood 
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(L: green), and the pseudo-likelihood (H: yellow; Huggins 
et al. 2010) when parameters values are N = 1000, and 
<p = 0.5 for p = 0.2 (top), and p = 0.5 (middle), p = 0.8 
(bottom). The long black horizontal lines show the 
expected population size at time 

Figure S7. Boxplots of abundance estimates N; for each 
sample time (k = 7) of 100 simulated datasets, for the 
Crosbie-Manley-Arnason-Schwarz (C: red), the likelihood 
(L: green), and the pseudo-likelihood (H: yellow; Huggins 
et al. 2010) when parameters values are N = 1000, and 
4> = 0.8 for p = 0.2 (top), and p = 0.5 (middle), p = 0.8 
(bottom). The long black horizontal lines show the 
expected population size at time 

Figure S8. Boxplots of estimated standard errors for the 
abundance estimates (SE(Nj)) for each sample time 
(k = 7) of 100 simulated datasets, for the Crosbie-Man- 
ley-Arnason-Schwarz (C: red), the likelihood (L: green), 
and the pseudo-likelihood (H: yellow; Huggins et al. 
2010) when parameters values are N = 200, and (p = 0.2 
for p = 0.2 (top), and p = 0.5 (middle), p = 0.8 (bottom). 
Estimates from simulations that produced a singular Hes- 
sian were removed. 

Figure S9. Boxplots of estimated standard errors for the 
abundance estimates (SE(Nj)) for each sample time 
(k = 7) of 100 simulated datasets, for the Crosbie-Man- 
ley-Arnason-Schwarz (C: red), the likelihood (L: green), 
and the pseudo-likelihood (H: yellow; Huggins et al. 
2010) when parameters values are N = 200, and <p = 0.5 
for p = 0.2 (top), and p = 0.5 (middle), p = 0.8 (bottom). 
Estimates from simulations that produced a singular Hes- 
sian were removed. 

Figure S10. Boxplots of estimated standard errors for the 
abundance estimates (SE(Nj)) for each sample time 
(k = 7) of 100 simulated datasets, for the Crosbie-Man- 



ley-Arnason-Schwarz (C: red), the likelihood (L: green), 
and the pseudo-likelihood (H: yellow; Huggins et al. 
2010) when parameters values are N = 200, and <p = 0.8 
for p = 0.2 (top), and p = 0.5 (middle), p = 0.8 (bottom). 
Estimates from simulations that produced a singular Hes- 
sian were removed. 

Figure Sll. Boxplots of estimated standard errors for the 
abundance estimates (SE(Nj)) for each sample time 
(k = 7) of 100 simulated datasets, for the Crosbie-Man- 
ley-Arnason-Schwarz (C: red), the likelihood (L: green), 
and the pseudo-likelihood (H: yellow; Huggins et al. 
2010) when parameters values are N = 1000, and <p = 0.2 
for p = 0.2 (top), and p = 0.5 (middle), p = 0.8 (bottom). 
Estimates from simulations that produced a singular Hes- 
sian were removed. 

Figure SI 2. Boxplots of estimated standard errors for the 
abundance estimates (SE(Nj)) for each sample time 
(k = 7) of 100 simulated datasets, for the Crosbie-Man- 
ley-Arnason-Schwarz (C: red), the likelihood (L: green), 
and the pseudo-likelihood (H: yellow; Huggins et al. 
2010) when parameters values are N = 1000, and <p = 0.5 
for p = 0.2 (top), and p = 0.5 (middle), p = 0.8 (bottom). 
Estimates from simulations that produced a singular Hes- 
sian were removed. 

Figure S13. Boxplots of estimated standard errors for the 
abundance estimates (SE(Nj)) for each sample time 
(k = 7) of 100 simulated datasets, for the Crosbie-Man- 
ley-Arnason-Schwarz (C: red), the likelihood (L: green), 
and the pseudo-likelihood (H: yellow; Huggins et al. 
2010) when parameters values are N = 1000, and <p = 0.8 
for p = 0.2 (top), and p = 0.5 (middle), p = 0.8 (bottom). 
Estimates from simulations that produced a singular Hes- 
sian were removed. 
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